---
title: "MTPS replication"
output:
  html_document: default
  pdf_document: default
date: "2022-11-14"
---
```{r, setup}

knitr::opts_knit$set(root.dir = "C:/Users/edmun/Documents/2023-2026 - DPhil Politics/2 - OTHER PAPERS/2 - UNDER REVIEW/UR - Twins and trust")
setwd("C:/Users/edmun/Documents/2023-2026 - DPhil Politics/2 - OTHER PAPERS/2 - UNDER REVIEW/UR - Twins and trust")
```

```{r Packages, echo = FALSE} 
library(OpenMx)
library(umx) 
library(haven)
library(tidyverse)
library(stargazer)
library(eeptools)
library(sjPlot)
library(ggplot2)
library(broom)
library(RColorBrewer)
library(psych)
library(ggh4x)
library(ggpubr)
set.seed(1234)
```

```{r Pre-Processing}

#Loading data
MTPS <- data.frame(read_spss("MTPS.sav"))


#Removing missing values and incomplete pairs
MTPS[MTPS == 999] <- NA
MTPS[MTPS == 99999999] <- NA
MTPS <- subset(MTPS, complete_pair == 1)

#Converting variables to factors
#MTPS[ , c(6:263, 268:271)] <- lapply(MTPS[ , c(6:263, 268:271)], as.factor)
MTPS$rsex <- as.numeric(MTPS$rsex)
MTPS$BDATE <- as.Date(MTPS$BDATE)
a = as.Date(c('2024-10-04'))
MTPS$AGE <- floor(age_calc(as.Date(MTPS$BDATE), enddate = a, units = "years", precise = TRUE)) - 14

#Turnout

MTPS$turnout04 <- as.factor(MTPS$Qq39)
MTPS <- MTPS %>% mutate(turnout04 = 
                          scale(dplyr::recode(turnout04, 
                                 "1" = 1,
                                 "2" = 0, 
                                 "3" = 0)))

MTPS$lifeturnout <- as.numeric(scale(6-MTPS$Qq41))

MTPS$overallturnout = as.numeric(scale(MTPS$turnout04 + MTPS$lifeturnout))

#Behavioral engagement

MTPS$rally <- as.factor(MTPS$Qq43_A_1)
MTPS <- MTPS %>% mutate(rally = 
                          as.numeric(scale(dplyr::recode(rally, 
                                 "1" = 1,
                                 "2" = 0))))

MTPS$workpol <- as.factor(MTPS$Qq43_A_2)
MTPS <- MTPS %>% mutate(workpol = 
                          as.numeric(scale(dplyr::recode(workpol, 
                                 "1" = 1,
                                 "2" = 0))))

MTPS$donate <- as.factor(MTPS$Qq43_A_3)
MTPS <- MTPS %>% mutate(donate = 
                          as.numeric(scale(dplyr::recode(donate, 
                                 "1" = 1,
                                 "2" = 0))))

MTPS$writepol <- as.factor(MTPS$Qq43_A_5)
MTPS <- MTPS %>% mutate(writepol = 
                          as.numeric(scale(dplyr::recode(writepol, 
                                 "1" = 1,
                                 "2" = 0))))

MTPS$behavioral <- as.numeric(scale(MTPS$donate + MTPS$writepol + MTPS$workpol + MTPS$rally))

#Political trust

MTPS$Qq27 = as.numeric(MTPS$Qq27)

MTPS$trustgov <- ((3-MTPS$Qq28))/2

MTPS$biginterest <- (MTPS$Qq29-1)

MTPS$trustpol <- (MTPS$Qq31-1)/2

MTPS$wastetax <- (MTPS$Qq30-1)/2

MTPS$trustgov_scaled = scale(MTPS$trustgov)

MTPS$biginterest_scaled = scale(MTPS$biginterest)

MTPS$wastetax_scaled = scale(MTPS$wastetax)

MTPS$trustpol_scaled = scale(MTPS$trustpol)

MTPS$trust <- as.numeric(MTPS$trustgov + MTPS$biginterest + MTPS$trustpol + MTPS$wastetax)/4
MTPS$trust_scaled = as.numeric(scale(MTPS$trust))

#Equal environment variables

MTPS$bedroom = as.numeric(MTPS$Qq73_A_1)
MTPS$dressalike = as.numeric(MTPS$Qq73_A_3)
MTPS$sameclass = as.numeric(MTPS$Qq73_A_4)
MTPS$samefriends = as.numeric(MTPS$Qq73_A_2)

#Creating twin ID variable
MTPS$TWINID = MTPS$COMPID %% 100
MTPS$TWINID = as.factor(MTPS$TWINID)
MTPS <- MTPS %>% mutate(TWINID = 
                          dplyr::recode(TWINID, 
                                 "0" = 1,
                                 "1" = 2))

#Stacking data by twin pairs
MTPS_desc = MTPS
MTPS <- umx_long2wide(MTPS, famID = "FAMID", twinID = "TWINID", zygosity = "ZYG")

#Creating MZ and DZ datasets
MTPS_MZ <- subset(MTPS, ZYG == 1)
MTPS_DZ <- subset(MTPS, ZYG == 2)

```

```{r Appendix 2: Descriptive statistics}

MTPS_desc %>% 
  select(AGE, ZYG, rsex, trust, trust_scaled, trustgov, wastetax, biginterest, trustpol, bedroom, dressalike, sameclass, samefriends, Qq27) %>% 
  filter(ZYG ==1) %>% 
  as.data.frame() %>% 
  rename(interpersonal = Qq27) %>% 
  select(!ZYG) %>% 
  select(!trust_scaled) %>% 
  stargazer(iqr = T, median = T, out = "MTPS_heritability_summary_MZ.htm")

MTPS_desc %>% 
  select(AGE, ZYG, rsex, trust, trust_scaled, trustgov, wastetax, biginterest, trustpol, bedroom, dressalike, sameclass, samefriends, Qq27) %>% 
  filter(ZYG == 2) %>% 
  as.data.frame() %>% 
  select(!ZYG) %>% 
  select(!trust_scaled) %>% 
  rename(interpersonal = Qq27) %>% 
  stargazer(iqr = T, median = T, out = "MTPS_heritability_summary_DZ.htm")

```

```{r Appendix 2: Political Trust ACE}

#Political trust
ACE <- umxACE("ACE", selDVs = "trust_scaled", selCovs = c("rsex", "AGE"), sep = "_T", mzData = MTPS_MZ, dzData = MTPS_DZ, optimizer = "SLSQP")
set.seed(1234)
ACE <- umxConfint(ACE, run = TRUE, level = .95, 
                  parm = c("top.A_std[1,1]", "top.C_std[1,1]", "top.E_std[1,1]"), 
                  wipeExistingRequests = F,
                  optimizer = "SLSQP")
summary(ACE, verbose = TRUE)

ACE1 <- umxACE("ACE1", selDVs = "trust_scaled", selCovs = c("rsex", "AGE", "sameclass", "samefriends", "bedroom", "dressalike"), sep = "_T", mzData = MTPS_MZ, dzData = MTPS_DZ, optimizer = "SLSQP")
set.seed(1234)
ACE1 <- umxConfint(ACE1, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.C_std[1,1]", "top.E_std[1,1]"), wipeExistingRequests = F, optimizer = "SLSQP")
summary(ACE1, verbose = TRUE)


```

```{r Appendix 2: Political trust AE and CE}
AE <- umxACE("AE", selDVs = "trust_scaled", selCovs = c("rsex", "AGE"), sep = "_T", mzData = MTPS_MZ, dzData = MTPS_DZ, optimizer = "SLSQP")
set.seed(1234)
AE = umxModify(AE, update = c("c_r1c1"), comparison = TRUE)
AE <- umxConfint(AE, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.E_std[1,1]"), wipeExistingRequests = F, optimizer = "SLSQP")
summary(AE, verbose = TRUE)

CE <- umxACE("CE", selDVs = "trust_scaled", selCovs = c("rsex", "AGE"), sep = "_T", mzData = MTPS_MZ, dzData = MTPS_DZ, optimizer = "SLSQP")
set.seed(1234)
CE = umxModify(CE, update = c("a_r1c1"), comparison = TRUE)
CE <- umxConfint(CE, run = TRUE, level = .95, parm = c("top.C_std[1,1]", "top.E_std[1,1]"), wipeExistingRequests = F, optimizer = "SLSQP")
summary(CE, verbose = TRUE)
```

```{r Appendix 2: Interpersonal trust heritability models}

ACET <- umxACE("ACE", selDVs = "Qq27", selCovs = c("rsex", "AGE"), sep = "_T", mzData = MTPS_MZ, dzData = MTPS_DZ, tryHard = "ordinal", optimizer = "SLSQP")
ACEt <- umxConfint(ACET, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.C_std[1,1]", "top.E_std[1,1]"), wipeExistingRequests = F, optimizer = "SLSQP")
summary(ACEt, verbose = TRUE)

ACEc <- umxACE("ACE", selDVs = "Qq27", selCovs = c("rsex", "AGE", "bedroom", "sameclass", "samefriends", "dressalike"), sep = "_T", mzData = MTPS_MZ, dzData = MTPS_DZ, optimizer = "SLSQP")
set.seed(1234)
ACEc <- umxConfint(ACEc, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.C_std[1,1]", "top.E_std[1,1]"), wipeExistingRequests = F, optimizer = "SLSQP")
summary(ACEc, verbose = TRUE)

AET = umxModify(ACET, update = c("c_r1c1"), comparison = TRUE)
AET <- umxConfint(AET, run = TRUE, level = .95, wipeExistingRequests = F, parm = c("top.A_std[1,1]", "top.E_std[1,1]"), optimizer = "SLSQP")
summary(AET, verbose = TRUE)

CET= umxModify(ACET, update = c("a_r1c1"), comparison = TRUE)
CET <- umxConfint(CET, run = TRUE, level = .95, wipeExistingRequests = F, parm = c("top.C_std[1,1]", "top.E_std[1,1]"), optimizer = "SLSQP")
summary(CET, verbose = TRUE)

umxCompare(ACET, comparison = c(AET, CET))

```

```{r Figure 2: Headline univariate models for all datasets}
univariate_results = read.csv("univariate_results.csv")

headline_figures = univariate_results %>% 
  filter(Model == "Political trust scale")

figure_All = headline_figures %>% 
  ggplot(
    aes(x = Parameter, y = Estimate) 
  ) + 
  geom_col(alpha = 0.3, 
           width = 0.5) + 
  geom_errorbar(aes(
    x = Parameter, 
    ymin = Lower, 
    ymax = Upper, 
  ), lty = 3, 
  width = 0.5) +
  ylim(0, 1) +
  theme_classic() +
  labs(x = "Parameter", y = "Standardised estimate") + 
  theme(text = element_text(size = 12, family = "serif")) + 
  facet_wrap(Dataset~., nrow = 1) + 
  scale_x_discrete(labels = c("Additive genetic", "Common environ.", "Unique environ."))

figure_All
```

```{r Figure 3: Interpersonal + institutional trust univariate models for all datasets}
univariate_results = read.csv("MTPS_univariate.csv")

headline_figures = univariate_results %>% 
  filter(Model == "Political trust scale" | Model == "Interpersonal trust" | Model == "Behavioral trust")

figure_All_interpersonal_US = headline_figures %>% 
  filter(Dataset == "United States (MN)") %>% 
  ggplot(
    aes(x = Parameter, y = Estimate) 
  ) + 
  geom_col(alpha = 0.3, 
           width = 0.5) + 
  geom_errorbar(aes(
    x = Parameter, 
    ymin = Lower, 
    ymax = Upper, 
  ), lty = 3, 
  width = 0.5) +
  ylim(0, 1) +
  theme_classic() +
  labs(x = "Parameter", y = "Standardised estimate") + 
  theme(text = element_text(size = 12, family = "serif")) + 
  facet_nested_wrap(Dataset~Model, nrow = 1) + 
  scale_x_discrete(labels = c("Additive Genetic.", "Common Environ.", "Unique Environ."))

figure_All_interpersonal_AUS = headline_figures %>% 
  filter(Dataset == "Australia") %>% 
  ggplot(
    aes(x = Parameter, y = Estimate) 
  ) + 
  geom_col(alpha = 0.3, 
           width = 0.5) + 
  geom_errorbar(aes(
    x = Parameter, 
    ymin = Lower, 
    ymax = Upper, 
  ), lty = 3, 
  width = 0.5) +
  ylim(0, 1) +
  theme_classic() +
  labs(x = "Parameter", y = "Standardised estimate") + 
  theme(text = element_text(size = 12, family = "serif")) + 
  facet_nested_wrap(Dataset~Model, nrow = 1) + 
  scale_x_discrete(labels = c("Add. Gen.", "Com. Env.", "Uniq. Env."))

figure_All_interpersonal_SE = headline_figures %>% 
  filter(Dataset == "Sweden") %>% 
  ggplot(
    aes(x = Parameter, y = Estimate) 
  ) + 
  geom_col(alpha = 0.3, 
           width = 0.5) + 
  geom_errorbar(aes(
    x = Parameter, 
    ymin = Lower, 
    ymax = Upper, 
  ), lty = 3, 
  width = 0.5) +
  ylim(0, 1)+
  theme_classic() +
  labs(x = "Parameter", y = "Standardised estimate") + 
  theme(text = element_text(size = 12, family = "serif")) + 
  facet_nested_wrap(Dataset~Model, nrow = 1) + 
  scale_x_discrete(labels = c("Additive Genetic.", "Common Environ.", "Unique Environ."))

figure_All_interpersonal = ggarrange(figure_All_interpersonal_AUS, figure_All_interpersonal_US, figure_All_interpersonal_SE, common.legend = T, legend = "right", ncol =1)

figure_All_interpersonal

``` 